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1.  Introduction 

Coastal  ocean  processes  are  complicated  and  they  happen  as  various  phenomena  that  span  a 
vast  range  of  spatial  and  temporal  scales.  For  instance,  general  circulations  of  oceans  occur 
at  global  scales  (Wunsch  &  Ferrari,  2004).  Tropical  waves  that  eventually  impact  coastal 
waters  propagate  with  wavelengths  of  one  thousand  kilometers  and  periods  of  one  month 
(Legeckis  et  al.,  1983).  Langmuir  cells,  which  are  pairs  of  vortices  hanging  below  water 
surfaces,  have  spatial  sizes  ranging  from  one  to  hundreds  of  meters  (Weller  et  al.,  1985). 
Scour  near  coastal  structures  is  under  influence  of  large-scale  current  processes  but  occurs  in 
a  relatively  small  size  (Sumer  &  Whitehouse,  2001).  Flows  around  enormous  numbers  of 
swimming  microorganisms  can  occur  at  scales  of  micrometers  (Pedley,  1992).  Here,  the 
scales  are  either  observation  scales  such  as  characteristic  length  and  time  or  process  scales 
such  as  those  in  wavelet  analysis  (Chui,  1992;  Kumar  and  Foufoula-Georgiou,  1997). 

Since  a  few  decades  ago,  various  geophysical  fluid  dynamics  (GFD)  models  have  been 
developed  for  individual  coastal  ocean  phenomena  at  specific  scales.  The  Princeton  Ocean 
Model  (POM),  Finite-Volume  Coastal  Ocean  Model  (FVCOM),  and  HYbrid  Coordinate 
Ocean  Model  (HYCOM)  were  developed  to  predict  current  velocity,  sea  level,  salinity,  and 
temperature  at  regional  scales  (Blumberg  and  Mellor,  1987;  Chen  et  al.,  2003;  Halliwell, 
2004).  The  WAVE  WATCH  and  SWAN  (Simulating  Waves  Nearshore)  models  were 
designed  to  simulate  surface  wave  propagation  at  global  to  coastal  scales  (Tolman,  1991; 
Booij  et  al.,  1999).  Models  have  also  been  proposed  to  predict  sediment  transport  and  seabed 
morphology  for  near-coastal  regions  (e.g.,  Tonnon  et  al.,  2007;  Papanicolaou  et  al.,  2008).  In 
recent  years,  computational  fluid  dynamics  (CFD),  which  can  accurately  model  small-scale 
and  detailed  flow  structures,  has  been  applied  to  coastal  engineering  flows  (e.g..  Young  et 
al.,  2001).  In  view  of  the  multiscale  and  multiphysics  nature  of  coastal  ocean  processes,  there 
is  a  great  challenge  to  simulate  them  accurately  and,  until  now,  the  efforts  using  numerical 
simulation  have  been  successful  merely  for  individual  phenomena  and  scales.  The  challenge 
comes  from  model  restrictions,  numerical  techniques,  and  computer  capabilities.  For 
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instance,  a  deep-ocean  model  has  difficulty  in  dealing  with  the  vertical  mesh  at  sudden 
bathymetry  changes  as  well  as  the  smaller  scales  of  nearshore  flows  (Song  &  Hou,  2006; 
Heimusund  &  Berntsen,  2004).  Limitations  such  as  hydrostatic  assumptions  and/or  two- 
dimensionality  of  GFD  models  are  inherent  restrictions  that  prohibit  accurate  simulations  of 
many  important  phenomena  such  as  vertical  motions  of  Langmuir  cells.  It  is  frequently 
reported  that  coastal  models  become  unstable  at  small  time  steps  and  grid  spacing  (e.g., 
Heimusund  &  Berntsen,  2004;  Keen  et  al.,  2003),  which  is  not  a  surprise  since  they  are 
designed  for  large-scale  flow  phenomena.  It  is  anticipated  that  the  solutions  of  the  models 
do  not  necessarily  converge  to  those  of  Navier-Stokes  equations  as  grid  spacing  tends  to 
zero.  Therefore,  it  is  not  realistic  to  achieve  accurate  simulation  of  coastal  phenomena, 
especially  the  small-scale  processes,  by  just  using  existing  coastal  ocean  models  and  simply 
reducing  time  steps  and  grid  spacing.  Although  in  principle  CFD  approaches  have  fewer 
limitations  and  can  capture  flow  phenomena  at  much  broader  spectra  of  scales,  they  are 
computationally  expensive  and  not  applicable  in  simulating  a  complete  actual  coastal  ocean 
flow.  It  should  be  noted  that,  although  both  GFD  and  CFD  are  based  on  the  Navier-Stokes 
equations,  they  are  different  approaches  with  respect  to  numerical  technique,  turbulence 
closure,  and  parameterization  for  small  scales. 

It  is  now  becoming  a  trend  in  prediction  of  coastal  ocean  flows  to  adapt  to 
multiphysics/ multiscale  approaches  (Fringer  et  al.,  2006).  Although  computer  modeling  has 
reached  the  point  where  the  simulation  of  individual  flow  phenomena  over  relatively 
narrow  ranges  of  scales  has  become  mature,  a  single,  comprehensive  model  capable  of 
dealing  with  multiphysics/ multiscale  problems  is  unlikely  to  be  available  in  the  near  future. 
The  hybrid  method  (HM)  couples  different  models  to  each  other,  and  the  domain 
decomposition  method  (DDM)  divides  a  flow  domain  into  many  subdomains,  each  of  which 
is  assigned  to  an  individual  model.  Combining  HM  and  DDM  is  one  of  the  most  promising 
currently  available  techniques  to  bridge  the  scales  and  overcome  difficulties  in 
multiphysics/ multiscale  modeling  (Benek  et  al.,  1983;  Harten,  1993;  Dolbow  et  al.,  2004). 
Since  erosion  and  transport  of  erodible  seabed  sediment  is  coupled  to  various 
hydrodynamic  forces,  it  is  imperative  to  analyze  them  at  different  scales  in  correspondence 
to  multiscale  hydrodynamic  processes  (Chiew,  1991).  It  has  been  common  to  assume  that 
seabed  scour  is  a  local  process,  which  occurs  within  a  few  tens  of  diameters  of  a  structure 
resting  on  the  seafloor  (Zang  et  al.,  2009).  This  is  true  for  event-scale  erosion  but  it  is 
unlikely  for  systematic  or  catastrophic  scour  and/  or  burial  processes  that  operate  at  months 
or  yearly  time  scales.  Experimental  results  and  parametric  methods  for  scour  around 
structures  on  the  seafloor  have  been  supplemented  by  numerical  models  that  focus  on  the 
finest  scales  around  the  obstacles  (Alam  &  Cheng,  2010).  These  high-resolution  models 
typically  focus  on  steady  flow  but  recent  laboratory  studies  have  examined  wave-induced 
scour  (Xu  et  al.,  2009)  and  numerical  methods  have  been  applied  to  shoaling  waves 
(Myrhaug  et  al.,  2008).  These  studies  use  simplified  hydrodynamic  forcing  because  of  the 
disparity  in  scales  between  the  external  (ocean  currents)  and  internal  flows  (around 
structures). 

CFD  approaches  for  the  fluid  flow  around  a  structure  are  not  easily  implemented  for 
sediments  because  of  the  interaction  of  turbulence  with  discrete  particles.  Thus,  scour 
models  have  traditionally  used  parametric  approaches  (Myrhaug  et  al.,  2008)  and  discrete 
particle  models  (Zamankhan  &  Doolatshahi,  2008).  Lattice  Boltzmann  methods  are  being 
explored  as  well  for  simple  geometries  (Alam  &  Cheng,  2010).  These  approaches  depend  on 
a  given  flow  field  that  is  typically  simplified  from  the  ocean  environment.  The  constitutive 
equations  for  these  models  are  derived  from  observations  in  flumes  and  the  seabed,  and 
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thus  are  specifically  formulated  for  near-field  flow.  The  other  extreme  is  to  use  a  numerical 
model  of  regional  flow  to  calculate  a  large-scale  mean  scour  that  can  be  used  as  an  indicator 
of  variations  in  potential  scour  within  a  selected  region  (Keen  &  Glenn,  2002).  These  models 
tend  to  use  the  same  macroscopic  approach  represented  by  the  Navier-Stokes  equations. 
There  is  a  lack  of  knowledge  of  the  impact  of  these  external  hydrodynamic  processes  on  the 
local  scour  and  burial  problem.  In  order  to  use  localized  scour  models  to  investigate  the 
impact  of  external  factors  on  seabed  scour,  it  is  necessary  to  implement  them  on  much  larger 
domains.  This  approach  is  problematic  with  current  computational  resources.  The  problem 
of  applying  a  macroscopic  model  (i.e.,  those  based  on  Navier-Stokes  equations)  for  the 
external  problem  and  a  scour  model  for  the  local  scour  is  prohibitive  because  of 
fundamental  differences  in  their  numerical  formulations.  This  problem  can  be  more  easily 
examined  using  macroscopic  approaches,  which  have  common  numerical  formulations  and 
parameterizations. 

This  chapter  summarizes  our  recent  work  in  modeling  of  multiscale  and  multiphysics 
hydrodynamics  phenomena  using  HM  and  DDM.  We  also  discuss  related  sediment  transport, 
with  emphasis  on  localized  scour  and  erosion  processes.  First,  a  hybrid  approach  that 
couples  the  FVCOM  and  a  CFD  model  is  described,  and  results  of  multiscale  simulation  for 
an  effluent  thermal  discharge  from  a  diffuser  at  the  ocean  bottom  is  presented.  Second,  an 
analysis  is  made  of  the  effects  of  local-scale  hydrodynamics  on  sediment  transport  around 
the  diffuser.  Third,  as  a  multiphysics  modeling  of  interaction  between  different  phenomena, 
simulation  of  flow  over  sand  dunes  under  action  of  surface  wind  is  presented  to  illustrate 
the  interaction  between  surface  waves,  currents,  and  morphology.  These  examples 
demonstrate  the  multiscale/ multiphysics  methodology  as  applied  to  problems  that  cannot 
be  simulated  as  either  local  or  external  phenomena.  They  also  indicate  that  multiscale  and 
multiphysics  simulations  of  hydrodynamics  are  more  advanced  than  conventional  modeling 
because  of  the  complex  interaction  between  the  flow  and  discrete  particles  in  the  ocean. 

2.  Hybrid  CFD  and  FVCOM  for  simulation  of  thermal  effluent  into  coastal 
flows 

A  well-tested  CFD  model  is  employed  in  this  work  (Lin  &  Sotiropoulos,  1997;  Tang  et  al., 
2003;  Tang  et  al.,  2008).  The  governing  equations  for  hydrodynamic  processes  of  the  CFD 
model  are  the  three-dimensional  (3D)  continuity  and  Navier-Stokes  equations  that,  in 
general  curvilinear  coordinates,  are  expressed  as  follows: 

r^+J^j(Fk-Fvk)  +  I  =  0,  (i) 

dt  d%k 

where 

T  =  diag(  0, 1, 1, 1),  Q  =  ( p,u ,  v,  w)T , 

Fk  =-(Uk,uUk  +  pg,vllk  +  p£,wUk  +  pg  )T , 
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where:  t  is  time;  /  =|  dx{  /  d%]  |  is  the  Jacobian  of  the  geometric  transformation  from 
Cartesian  coordinates,  x;  (i  =  1,  2,  3  for  x,  y,  z  axes,  respectively)  to  curvilinear  coordinates 
( j  =  1,  2,  3);  p  is  the  static  pressure  divided  by  the  density;  u,  v,  and  w  are  the  velocities  in 
x,  y,  and  z  directions,  respectively;  T  is  the  temperature.  Furthermore,  Uk  =  u£xx^  are  the 
contravariant  velocities  in  £ k  directions  (Note  that  m  are,  respectively,  u,  v,  and  w,  g  are 
the  metrics  of  the  geometric  transformation,  and  repeated  indices  imply  summation);  Re  is 
the  Reynolds  number;  vtis  turbulence  eddy  viscosity;  gl]  ( gl]  =  %lx^  )  is  the  contravariant 
metric  tensor;  Fr  is  Froude  number;  e  is  the  unit  in  the  gravity  direction.  The  governing 
equation  for  heat  transfer  reads  as 

1ST _ 8  \l(  1  Vt_ 

]  dt  8^k  | /  ( Pr Re  +  Prf  y  } 

where:  Pr  is  the  molecular  Prandtl  number;  Prt  is  the  turbulent  Prandtl  number;  l  =  1,  2,  3. 
The  standard  mixing  length  model  is  used  in  this  work  (Mason,  1989).  The  governing 
equations  are  discretized  using  a  second-order  accurate,  implicit,  finite-volume  method  on 
non-staggered  grids,  and  they  are  solved  using  a  dual  time-stepping  artificial  compressibility 
method.  A  fourth-difference  artificial  dissipation  method  and  V-cycle  multigrid  method  are 
used.  A  domain  decomposition  method  in  conjunction  with  the  Schwartz  alternative 
iteration  is  employed  to  deal  with  complicated  geometry.  In  order  to  achieve  seamless 
transition  of  solutions  between  subdomains,  an  effective  mass  conservation  algorithm  is 
proposed.  The  CFD  model  has  been  tested  and  applied  in  various  problems  from  academe 
as  well  as  industry,  such  as  vortex  breakdown  and  flow  past  bridge  piers  (Sotiropoulos  and 
Ventikos,  1998;  Ge  and  Sotiropoulos,  2005).  For  details  of  the  CFD  model,  readers  are 
referred  to  Lin  &  Sotiropoulos  (1997),  Tang  et  al.  (2003;  2008),  and  Xu  &  Sun  (2009). 

In  the  FVCOM,  the  governing  equations  are  the  continuity  and  momentum  equations.  The 
governing  equations  for  the  external  mode  in  the  model  are  vertically  averaged  two- 
dimensional  (2D)  continuity  and  momentum  equations.  The  2D  continuity  equation  is  (Chen 
et  al.,  2006): 


lk  dT 
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where:  rj  is  the  water  surface  elevation;  H  is  water  depth;  and  Ui  are  depth-average  current 
velocities  in  xi  directions.  The  momentum  equations  are: 
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where:  the  overbar  denotes  the  vertical  integration;  g  is  the  gravity  constant;  po  is  the 
reference  density  of  seawater;  a  is  vertical  coordinate;/  is  the  Coriolis  parameter;  rsx,  are 
the  surface  friction  stresses,  and  rhx,  are  the  bottom  friction  stresses  in  Xi  directions;  D  is  the 
mean  water  depth;  and  Am  is  the  horizontal  eddy  viscosity  (Smagorinsky,  1963).  The 
governing  equations  of  the  internal  mode  in  the  model  are  3D  continuity  and  momentum 
equations  with  x  and  y  as  horizontal  coordinates  and  a  as  vertical  coordinate: 
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where:  m  are  the  layer  horizontal  velocity  components;  co  is  the  vertical  velocity  in  a 
coordinate;  and  Km  is  the  vertical  eddy  viscosity,  which  is  determined  using  the  Mellor  and 
Yamada  level-2.5  turbulent  closure  (Mellor  &  Yamada,  1982;  Chen  et  al.,  2003).  In  the  FVCOM: 
i,  j,  1=1,  2  (i  *  j  ).  The  governing  equations  are  discretized  using  a  finite-volume  method  with 
triangle  meshes  on  horizontal  planes  and  layer  meshes  in  the  vertical  direction.  Second-order 
accurate  upwind  schemes  are  used  to  discretize  the  advection  terms,  and  Runge-Kutta 
methods  are  employed  to  march  in  time.  The  external  and  internal  modes  may  have  different 
time  steps.  Mainly  because  it  uses  unstructured  meshes,  the  FVCOM  is  becoming  popular  in 
coastal  ocean  modeling.  Details  for  the  FVCOM  can  be  found  in  Chen  et  al.  (2003;  2006). 

In  order  to  simulate  multiscale  and  multiphysics  coastal  hydrodynamic  processes,  we 
proposed  HM  and  DDM  approaches  (Tang  &  Wu,  2010;  Wu  &  Tang  2010).  In  particular,  the 
CFD  model  is  coupled  with  the  FVCOM;  the  CFD  model  is  employed  to  resolve  small-scale 
flow  phenomena,  and  the  FVCOM  is  used  to  model  background  circulation.  The  domains  of 
the  CFD  model  and  FVCOM  overlap  over  a  region  (Fig.  1).  As  a  coupling  strategy,  the  3D 
CFD  model  is  coupled  to  the  3D  internal  mode  of  the  FVCOM,  and  the  two  models 
exchange  solutions  for  the  velocity  distributions  at  grid  interfaces  between  them.  The 
strategy  is  based  on  the  assumption  that  the  obstacle  covered  by  the  CFD  model  only  alters 
local  velocity  distribution  but  does  not  affect  horizontal  average  velocities  and  water  surface 
elevation,  which  are  determined  by  the  external  mode.  This  assumption  is  consistent  with 
the  assumption  in  the  FVCOM  (Chen  et  al.,  2003).  Since  Chimera  grids  (grids  overlapping 
arbitrarily  with  each  other)  provide  the  best  possible  flexibility  in  connecting  different 
models,  we  employed  them  to  connect  the  CFD  model  and  the  FVCOM,  as  shown  in  Fig.  2. 
In  the  figure,  a-b-c  and  d-e-f  are  grid  interfaces  between  CFD  model  and  FVCOM,  and  linear 
interpolations  are  made  to  provide  solutions  for  velocities  at  the  nodes  and  elements  on  the 
interfaces.  The  details  of  the  coupling  methods  can  be  found  in  Tang  &  Wu  (2010)  and  Wu  & 
Tang  (2010). 


76 


Sediment  Transport 


wave 


FVCOM  FVCOM 


undersea  obstacle 


Fig.  1.  Schematic  representation  of  CFD  model  and  FVCOM  coupling. 


c 


Fig.  2.  Layout  of  grids  of  CFD  and  FVCOM.  Solid  lines  -  FVCOM  grid,  dash  lines  -  CFD 
grid,  a)  Horizontal  plane  grids,  b)  Vertical  direction  grids. 
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Fig.  3.  Region  of  study,  bathymetry,  FVCOM  mesh,  and  diffuser  location. 
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Consider  an  effluent  discharge  from  a  diffuser  in  the  setting  of  the  New  York  and  New 
Jersey  coastal  region  under  action  of  tides  (Fig.  3).  The  x  axis  is  in  the  east  direction,  and 
y  axis  is  in  the  north  direction.  The  diffuser  is  located  in  the  New  York  Bight.  The  diffuser 
consists  of  a  pipe  with  a  diameter  of  1.32  m  lying  on  the  ocean  bottom,  and  ten  ports  with 
diameters  of  0.175  m  (Fig.  4a).  All  of  the  ten  ports  have  an  angle  of  110°  with  respect  to  the 
x  axis  and  upward  angles  with  respect  to  the  x-y  plane  ranging  evenly  from  45°  to  18°  from 
the  first  port,  located  near  the  origin  of  the  plane,  to  the  last  port.  Hot  water  at  32.0  °C  is 
discharged  at  speed  3.92  ms4  from  the  ports  into  the  coastal  water  at  20.5  °C.  The  only 
driving  force  of  the  flow  is  tides;  wind  and  other  factors  are  ignored.  This  is  a  multiscale 
flow;  the  thermal  effluent  happens  at  scales  of  the  discharge  ports  and  the  pipe,  while  the 
unsteady  tides  occur  at  scales  of  the  tides. 

The  hybrid  approach  is  used  to  simulate  the  flow;  the  CFD  model  captures  the  flow  around 
the  effluent  from  the  diffuser,  and  the  FVCOM  describes  the  background  large-scale  flows. 
The  mesh  of  the  FVCOM  is  shown  in  Fig.  3,  and  those  of  the  CFD  model  for  the  diffuser  are 
shown  in  Fig.  4b.  The  CFD  model  and  modeling  of  such  thermal  discharges  have  been 
intensively  tested  and  calibrated  using  other  models  and  measurement  data  (Tang  et  al., 
2008).  The  coupling  approach  has  been  also  validated  (Tang  &  Wu,  2010;  Wu  &  Tang,  2010). 
The  computed  solutions  are  shown  in  Figs.  5,  6,  7,  and  8.  Fig.  5  presents  the  large-scale 
coastal  flows  at  ebb  and  flood  tides  obtained  with  the  FVCOM.  It  is  seen  from  the  figure  that 
there  are  many  total  velocity  patches,  in  red  and  blue,  at  scales  ranging  from  104  to  105  m. 
The  thermal  discharge,  located  near  x,y  =  0,  is  at  the  edges  of  a  relative  high  velocity  patch. 
No  solution  details  for  flows  around  the  diffuser  are  available  from  simulations  of  the 
FVCOM;  however,  the  CFD  model  provides  the  details.  For  example.  Fig.  6  presents  velocity 
distribution  at  a  plane  6  m  above  the  diffuser.  Due  to  the  presence  of  the  diffuser  and 
thermal  effluent,  the  flow  field  during  both  ebb  and  flood  tides  is  greatly  altered;  there  are 
several  velocity  patches  at  scales  of  10  m  and  larger,  together  with  a  low-velocity  region  at 
downstream  side  of  the  diffuser,  which  runs  from  x,y  -  0  to  the  northwest  direction  (Fig.  6). 
Especially,  there  is  a  large  vortex  right  behind  the  diffuser  in  case  of  flood  tides. 


Fig.  4.  (a)  Configurations  of  the  diffuser,  (b)  CFD  Meshes  for  the  diffuser 
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As  seen  in  Fig.  4a,  the  discharge  ports  point  to  the  northeastern  direction,  and  thus  the  hot 
water  jets  should  be  in  the  same  direction.  Fig.  7  indicates  that  during  ebb  tides,  the  jets  are 
initially  towards  northeast  but  shortly  later,  due  to  the  northwestern  ambient  tide  currents, 
they  turn  towards  the  northwestern  direction.  While  during  flood  tides,  the  jets  and  the 
currents  are  in  about  same  directions,  the  thermal  plume  runs  to  a  far  downstream  location. 
The  CFD  model  has  a  high  resolution  at  the  mouth  of  the  ports,  with  10  grid  nodes  across 
the  port  diameter  of  order  1  cm  resolution.  As  demonstrated  in  our  previous  studies  (Tang 
et  al.,  2008),  the  CFD  model  accurately  resolves  velocities  and  temperature  at  the  mouths  of 
all  ports  using  this  mesh  resolution. 
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Fig.  5.  Large-scale  surface  velocity  field  at  (a)  flood  tides  and  (b)  ebb  tides. 
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Fig.  6.  Local  velocity  on  a  plane  6  m  above  the  diffuser  at  (a)  flood  and  (b)  ebb  tide. 
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(a) 


Fig.  7.  3D  temperature  plume  at  (a)  flood  and  (b)  ebb  tide. 
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Fig.  8.  Cross-section  velocity  and  temperature  at  10  m  downstream  from  the  diffuser  in  the  x 
direction,  (a)  Flood  tide,  (b)  Ebb  tide. 

Figure  8  shows  the  temperature  and  velocity  field  at  a  vertical  cross-section,  which  is  at  the 
east  side  of  diffuser.  In  case  of  the  flood  tides,  the  currents  and  thermal  jets  from  all  ports  are 
in  opposite  directions  and  the  thermal  plume  cannot  reach  this  section  whereas,  for  ebb 
tides,  they  are  in  the  same  directions  and  the  plume  is  still  strong  on  the  section.  In  the  latter 
case,  the  jets  from  individual  ports  already  have  merged  with  each  other  and  no  individual 
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plumes  are  observed,  and  this  is  the  result  of  intensive  mixing  in  the  cross-section  (Fig.  8b). 
As  shown  in  Fig.  8,  there  are  fast  lateral  currents  along  the  cross-section.  Due  to  interaction 
of  the  effluent  flows  and  the  cross-section  currents,  there  are  vertical  vortices  next  to  the 
head  of  the  diffuser,  located  near  x,  y  =  0.  It  should  be  noted  that  all  these  small  structures  of 
the  flow  and  temperature  fields  are  only  captured  by  the  CFD  model,  and  they  are  not 
available  from  the  FVCOM. 

3.  Scour  around  a  diffuser  at  the  seabed 

This  section  discusses  erosion  and  scour  around  the  diffuser  based  on  the  modeling  results 
for  the  multiscale  hydrodynamic  processes  in  the  previous  section.  Including  potential  scour 
computations  on  all  of  these  FVCOM  and  CFD  grids  is  beyond  the  scope  of  the  current 
work.  Instead,  the  results  of  the  flow  field  at  1  mab  (meter  above  bed)  from  the  overlapping 
grids  are  averaged  for  1  m2  boxes.  The  area  of  interest  is  centered  on  the  diffuser  and 
represented  by  a  grid  of  40  x  70  cells  that  are  1  m  along  each  axis.  The  output  from  the 
different  grids  is  available  at  a  range  of  heights  above  the  seafloor;  we  have  chosen  1  m  as 
the  representative  height  for  the  shear  stress  computations.  All  of  the  grids  are  included  in 
the  average,  which  produces  fields  for  the  flood  and  the  ebb  tides. 

No  wave  data  are  used  in  these  hydrodynamic  simulations,  which  aim  only  to  show  the 
influence  of  tidal  flow.  The  sedimentation  computations  do  not  include  advection  and  no 
sediment  is  transported;  the  plots  thus  show  the  equilibrium  resuspension  depth  Hr  as 
defined  by  Keen  and  Glenn  (1998).  These  depths  are  a  good  indication  of  scour  at  shorter 
time  scales.  The  scour  depths  shown  in  the  plots  represent  the  total  that  would  occur  over 
an  approximately  1-hour  period  for  which  the  bottom  currents  would  be  representative. 

The  mean  currents  during  flood  tide  are  very  low  but  the  distribution  has  additional 
structure  not  seen  in  the  FVCOM  output  (Fig.  9a).  A  very  similar  pattern  is  predicted  during 
ebb  tide  (Fig.  9b).  This  area  is  likely  to  be  the  focus  for  erosion.  The  bottom  sediment  at  this 
general  location  is  expected  to  be  variable  because  of  its  proximity  to  the  Hudson  River 
Canyon  and  the  resulting  large  sediment  outflow  from  the  Hudson  River.  The  available 
sediment  samples  indicate  available  sediment  from  clay  to  sand  (USGS  2010).  We  are  using 
a  representative  sediment  distribution  of  20  classes  ranging  from  4  microns  to  2000  microns 
(clay  to  gravel).  The  uniform  sediment  is  represented  by  a  mean  of  100  microns  (very  fine 
sand)  with  a  standard  deviation  of  40  microns.  This  specific  sediment  population  includes 
55%  clay  and  44%  very  fine  sand,  with  <  1  %  fine  sand.  The  computed  critical  shear  stresses 
are  0.465  and  0.4388  Pa,  respectively  for  the  dominant  sizes. 

Detailed  simulations  of  erosion  would  require  unavailable  field  data  for  the  area.  We  are 
using  the  erosion  rate  formulation  from  Ariathurai  et  al  (1983): 


El  =  M(rb/  rC/i  -1)  (10) 

where:  M  is  the  entrainment  rate;  n  is  the  bottom  shear  stress;  rC/i  is  the  critical  shear  stress 
for  entrainment  of  size  class  i.  This  formula  is  useful  for  adjusting  entrainment  for  unknown 
areas.  The  value  of  M  used  in  these  simulations  is  0.001  kg  nr2  s4.  The  bottom  shear  stress  is 
estimated  from  the  common  formulation: 


Tb  =  A,Q(u2+v2)  (11) 

where:  pw  is  the  water  density,  which  is  set  tol025  kg  nv3,  and  Q  is  the  coefficient  of  drag, 
which  is  set  to  1.  The  large  drag  coefficient  is  necessary  because  the  currents  from  the  model 
were  too  low  to  exceed  rc  for  the  sediment  used.  Including  wave  effects  would  improve  this 
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calculation.  The  resulting  field  of  %  shows  the  effect  of  the  diffuser  (Fig.  10).  Note  that  the 
bed  stress  is  not  symmetrical.  Larger  shear  stresses  are  generated  during  the  ebb  tide  and 
the  pattern  is  slightly  different  from  the  flood.  The  persistent  high  shear  stress  is  caused  by 
the  increased  flow  around  the  diffuser.  This  is  partly  attributable  to  the  consistent  seaward 
flow  near  x  =  -8  m  during  both  flood  and  ebb  tides,  which  is  caused  by  a  diffuser  outflow  of 
~3.92  ms4.  The  diffuser  flow  is  offshore,  which  appears  to  reinforce  the  ebb  tide  and  make  it 
more  localized.  This  results  in  the  concentration  of  higher  stress  around  the  outlets. 


X  (m)  X  (m) 

(a)  (b) 


Fig.  9.  Averaged  currents  from  all  CFD  grids  at  1  mab.  (a)  Flood  tide;  maximum  velocity  is 
0.40223  m  s4.  (b)  Ebb  tide;  maximum  velocity  is  0.40676  m  s-1. 
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(a)  (b) 

Fig.  10.  Bottom  drag  calculated  from  Eq.  (11).  (a)  Flooding  tide;  maximum  shear  stress  is 
165.83  Pa.  (b)  Ebb  tide;  maximum  shear  stress  is  169.59  Pa. 
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The  pattern  of  potential  scour  for  the  flood  (Fig.  11a)  and  ebb  (Fig.  lib)  tides  is  reflective  of 
the  shear  stresses.  Note  the  greater  scour  seaward  of  the  northern  end  of  the  diffuser  during 
the  ebb  tide.  This  pattern  is  due  to  the  fact  that,  as  shown  in  Fig.  4a,  the  ports  near  the  end 
have  lower  angles  in  comparison  with  the  ports  at  the  other  end,  and  this  leads  to  a  stronger 
horizontal  velocity  and  more  intensive  shearing  at  the  end.  These  calculations  of  potential 
scour  demonstrate  the  importance  of  calculating  the  detailed  flow  around  the  obstruction. 
No  calculations  were  performed  for  the  FVCOM  flow  fields  because  it  provided  uniform 
near-bottom  flow  with  a  small  magnitude.  These  simulations  suggest  that  there  could  be 
substantial  scour  and  possible  damage  to  the  structure  if  not  considered  in  its  design.  The 
scour  depths  would  be  greatly  increased  by  the  impact  of  storm  waves  1-4  m  high,  which 
occur  frequently  during  the  fall  and  winter  in  the  Middle  Atlantic  Bight  (Keen  &  Glenn, 
1995). 
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Fig.  11.  Potential  scour  calculated  from  Eq.  (10).  (a)  Flood  tide;  maximum  erosion  is  370.35 
mm.  (b)  Ebb  tide;  maximum  erosion  is  378.76  mm. 


4.  Evolution  of  sand  dunes  under  action  of  surface  waves 

The  governing  equation  for  wave  action  is  given  in  a  conservation  form  as  (e.g.,  Mei,  1983): 

8_N  |  dC^N  |  sc^N  |  5CgN  S 
dt  dxl  da  dO  a  ’ 

where:  t  is  the  time;  l  =1,  2,  and  xi  correspond  to  x  and  y  directions;  a  is  the  frequency;  and 
6  is  the  angle  of  the  wave  propagation  direction;  N  is  the  wave  action;  CgXj  are  the  wave 
speed  in  xi  direction  in  the  physical  space  ( x ,  y);  Ca  and  Ce  are,  respectively,  wave  speed  in 
a  and  6  direction  in  the  spectrum  space  ( a ,  6 );  S  is  a  source  term  that  represents  the 
combined  effects  of  wind  and  other  processes.  In  the  current  study,  short  waves  are 
considered,  which  leads  to  a  simple  expression  for  CgX/  (Mei,  1983). 
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The  shallow  water  equations  consist  of  the  equation  of  mass  conservation  (e.g.,  Mei,  1983) 


DH  dHUi 
- + - - 

dt  dxx 


=  0, 


and  the  equations  of  momentum  conservation: 


(13) 
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where:  i,l  =1,  2;  vt  is  the  turbulence  eddy  viscosity;  p  is  the  density;  rj  is  the  water  surface 
elevation;  Sx,x  are  the  radiation  stresses  resulting  from  waves. 

The  seabed  morphology  is  controlled  by  the  Exner  equation  (e.g.,  Henderson,  1966;  Sleath, 
1984): 


8H^+dBU1(Ul  +  Lpl  Q 

dt  dxx  K 

Here  Hb  is  the  elevation  of  the  seabed,  and  B  and  q  are  constants.  Typically,  B  depends  on 
flow  velocity,  water  depth,  sediment  grain  sizes,  and  other  factors,  and  q  is  usually  in  the 
range  of  0.5  <  q  <  1.5  .  In  this  paper,  q  =  1  is  used  (Hudson  &  Sweby,  2005;  Kubatko  et  al., 
2006). 

The  interactions  among  waves,  currents,  and  morphology  can  be  seen  in  the  governing 
equations.  In  Eq.  (12),  the  wave  action  is  related  to  wave  speed  CgX;  ;  these  are  coupled  with 
the  current  through  the  velocity  Ui  in  Eq.  (12)  (Mei,  1983).  As  indicated  in  Eq.  (14),  the 
current  is  affected  by  the  wave  field  through  the  radiation  Sxx  ,  and  morphology  evolution 
through  bottom  elevation  Hb  (77  =  H  +  Hb).  Actually,  the  current  is  also  affected  by  wind 
through  the  bottom  stresses  rbx .  .  Morphology  is  directly  related  to  the  current  through  the 
velocity  field  as  shown  in  Eq.  (15).  Details  for  the  governing  equations  can  be  found  in  Tang 
etal.  (2009). 

Equations  (12)  through  (15)  comprise  a  coupled  non-homogeneous  system  of  conservation 
laws.  Each  component  in  the  system  reproduces  the  framework  of  well-known  models.  For 
instance,  the  wave  action  equation  (12)  is  employed  in  SWAN  (Booij  et  al.,  1999),  the  shallow 
water  equations  (13)  and  (14)  are  used  for  SHORECIRC  (Luettich  and  Westerink,  2004),  and 
the  morphology  evolution  equation  (15)  is  also  widely  employed  in  engineering  (e.g.,  Wu, 
2004;  Kubatko  et  al.,  2006).  In  order  to  solve  the  system,  an  extension  of  the  Lax-Friedrich 
scheme  (Lax,  1954)  is  applied  to  discretize  the  wave  action  equation  (12).  Second,  the 
MacCormack  scheme  (MacCormack,  1969)  is  employed  to  solve  the  shallow  water 
equations.  Finally,  the  fourth  Euler  scheme  and  central  difference  operator  are  used  to  solve 
the  morphology  equation.  The  coupled  system  and  the  code  have  been  validated  and 
calibrated  using  a  series  of  problems.  For  details  of  the  discretization,  validation,  and 
calibration,  the  readers  are  referred  to  Tang  &  et  al.  (2009). 

The  above  coupled  wave,  current,  and  morphology  system  is  applied  to  study  evolution  of 
sand  dunes  under  action  of  surface  waves  on  the  horizontal  plane.  The  initial  conditions  for 
the  wave  energy,  the  velocity  field,  and  the  bottom  shape  are,  respectively 


N= 0, 


(16a) 
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Ui=0,  U2=0,  H=2- Jexp{-0.01((x-0.5f)2  +y2)} ,  (16b) 

i= 1 

^  exp  j-0.01((x-  0.5m)2  +  y2)J  /  (16c) 

m=l 

where  n  is  the  number  of  sand  dunes.  In  this  study,  n  =  1,  2,  3.  The  upstream  boundary 
condition  is 

H  =  2m,  x  =  -15m.  (17) 

The  wind  effect  is  introduced  as  a  source  term  in  Eq.  (12)  at  the  upstream  end: 

S  -  0.15(1  -  tanh(20  +  x  -  O.Olf )) ,  (18) 

which  generates  surface  waves  starting  at  the  upstream  end  and  propagating  in  the  x 
direction.  In  the  simulation.  Ax,  Ay  =0.4  m.  Act  =  0.1  s4,  and  A 0  =0.76  radians.  Extrapolation 
is  used  to  update  solutions  at  the  boundaries.  Other  related  parameters  are  the  same  as 
those  presented  in  our  previous  study  on  a  single  sand  dune  case  (Tang  et  al.,  2009). 


(b)  (c) 

Fig.  12.  Evolution  of  sand  dunes  (seabed  elevation) 

As  seen  in  Eq.  (18),  the  driving  force  S  is  a  wave  propagating  in  the  x  direction.  The  driving 
force  generates  a  train  of  surface  waves  traveling  in  the  same  direction.  The  computed 
instantaneous  evolution  of  sand  dune  elevation,  surface  height,  and  wave  action,  are 
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presented  in  Figs.  12  through  14  for  three  sand  dunes  (n= 3).  As  shown  in  Fig.  12,  as  the  wind 
blows  from  left  to  the  right,  all  three  sand  dunes  move  in  the  x  direction,  and  gradually  they 
merge  with  each  other.  During  this  process,  each  sand  dune  is  changing  from  a  circle  into  a 
triangle,  a  steroidal,  and  then  a  strip,  which  are  the  three  typical  shapes  of  sand  dunes  in 
their  evolution.  At  the  same  time,  low-elevation  regions  are  forming  at  the  front,  lateral 
sides,  and  immediately  behind  the  three  sand  dunes.  In  Fig.  13,  it  is  seen  that  the  water 
surface  has  a  bump  in  front  and  a  dip  behind  each  of  the  three  dunes.  In  comparison  with 
sand  dune  configurations  in  Fig.  12,  it  is  known  from  Fig.  14  that  wave  action  reaches 
minimums  above  all  sand  dunes  and  maximums  right  behind  them.  From  Figs.  12  and  13  it 
is  seen  that,  during  the  evolution  of  the  sand  dunes,  both  the  wave  field  and  water  surface 
elevation  evolve  congruently  with  the  sand  dunes.  Figure  14  reveals  clear  traces  of  bed 
morphology,  which  indicates  its  strong  effect  on  the  wave  field. 
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Fig.  13.  Evolution  of  water  surface  elevation 

Figure  15  shows  the  evolution  of  the  heights  and  locations  of  three  sand  dunes  as  the  wave 
is  propagating  towards  the  right.  An  interesting  note  is  that  the  height  of  the  first  dune,  the 
leading  dune  that  is  most  upstream,  decreases  faster  at  the  beginning  (t  <  470  s)  in  comparison 
with  the  heights  of  the  other  two  dunes.  As  seen  in  Fig.  15a,  however,  at  about  t  =  470  s  the 
second  and  the  third  disappear  or  they  merge  with  the  first  one,  whereas  the  first  still  exists. 
After  the  disappearance  of  the  second  and  third  sand  dunes,  the  height  of  the  first  dune 
remains  about  the  same  for  a  while.  It  is  seen  in  Fig.  15b  that,  as  all  three  dunes  are  about  to 
merge  with  each  other  at  t  =  470  s,  the  first  sand  dune  moves  faster  than  the  other  two,  with 
a  large  dx/dt,  which  accelerates  the  merging  process. 
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Fig.  14.  Evolution  of  wave  action 


Fig.  15.  Evolution  of  3-dune,  a)  Heights  of  sand  dunes,  b)  Location  of  sand  dunes 

Comparison  of  Figs.  12  through  14  with  those  obtained  from  modeling  of  single  sand  dunes 
(not  shown)  suggests  that  the  evolution  of  a  sand  dune  in  the  case  of  multiple  sand  dunes  is 
consistent  with  a  single  sand  dune  in  certain  aspects,  such  as  the  shape  of  sand  dunes 
(e.g.,  Hudson  &  Sweby,  2005;  Tang  et  al.,  2009).  Nevertheless,  sand  dune  development  in 
situations  of  multiple  sand  dunes  indeed  behaves  differently.  In  Fig.  16,  the  simulated  results 
for  situations  of  1,  2,  and  3  sand  dunes  (n  =  1,  2,  and  3  in  Eq.  16)  are  presented,  which  shows 
the  heights  and  locations  of  the  leading  or  the  most  upstream  sand  dunes.  The  most  interesting 
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finding  in  this  figure  is  that,  in  case  of  multiple  sand  dunes,  the  heights  of  the  leading  sand 
dunes  decrease  slower  but  they  move  downstream  faster.  The  larger  the  number  of  the  sand 
dunes,  the  more  pronounced  this  trend.  It  is  also  seen  that  the  evolutions  of  heights  and 
locations  of  the  leading  sand  dunes  in  case  of  2  and  3  sand  dunes  are  not  much  different,  but 
they  are  obviously  distinct  from  that  of  an  isolated  sand  dune  (Fig.  16).  This  is  a  clear 
indication  of  the  interaction  between  individual  sand  dunes  and  their  resulting 
hydrodynamics. 


Fig.  16.  Evolutions  of  the  leading  sand  dunes  in  case  of  1-,  2-,  and  3-dune,  a)  Heights  of  sand 
dunes,  b)  Locations  of  sand  dunes 

5.  Concluding  remarks 

This  chapter  presents  the  needs  to  simulate  multiscale  and  multiphysics  coastal  ocean 
hydrodynamics  and  the  necessity  to  include  its  effects  in  estimation  of  seabed  sediment 
scour  and  morphology  evolution.  It  describes  the  multiscale  and  multiphysics  method  the 
authors  proposed  as  a  prominent  approach,  which  is  a  hybrid  method  in  conjunction  with 
the  domain  decomposition  method.  The  feasibility  and  potential  of  the  approach  in  resolving 
multiscale  and  multiphysics  processes  is  demonstrated  using  example  computations. 

The  potential  scour  for  the  seabed  diffuser  indicates  two  important  results:  (1)  the  additional 
detail  of  the  flow  field  computed  by  the  CFD  model  is  critical  in  capturing  the  areas  of  scour 
by  tides;  and  (2)  it  is  necessary  to  include  waves  for  both  mean  and  storm  conditions  to 
estimate  the  cumulative  potential  scour  around  the  diffuser.  The  local  scour  problem  can  be 
computed  using  a  high-resolution  model  but  this  would  not  permit  an  examination  of 
variations  within  the  seafloor  area  of  interest.  This  is  especially  true  for  waves,  which  are 
sensitive  to  slight  changes  in  water  depth.  The  modeling  of  sand  dunes  illustrates  the  strong 
interaction  between  surface  waves,  currents,  and  morphology.  The  simulation  presents 
interesting  features  of  sand  dunes  with  respect  to  their  heights  and  locations,  and  it  clearly 
indicates  that  evolution  of  multiple  sand  dunes  is  distinct  from  that  of  single  sand  dunes. 
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